/*
 * vim:syntax=c
 * A small module to implement missing C99 math capabilities required by numpy
 *
 * Please keep this independant of python ! Only basic types (npy_longdouble)
 * can be used, otherwise, pure C, without any use of Python facilities
 *
 * How to add a function to this section
 * -------------------------------------
 *
 * Say you want to add `foo`, these are the steps and the reasons for them.
 *
 * 1) Add foo to the appropriate list in the configuration system. The
 *    lists can be found in numpy/core/setup.py lines 63-105. Read the
 *    comments that come with them, they are very helpful.
 *
 * 2) The configuration system will define a macro HAVE_FOO if your function
 *    can be linked from the math library. The result can depend on the
 *    optimization flags as well as the compiler, so can't be known ahead of
 *    time. If the function can't be linked, then either it is absent, defined
 *    as a macro, or is an intrinsic (hardware) function.
 *
 *    i) Undefine any possible macros:
 *
 *    #ifdef foo
 *    #undef foo
 *    #endif
 *
 *    ii) Avoid as much as possible to declare any function here. Declaring
 *    functions is not portable: some platforms define some function inline
 *    with a non standard identifier, for example, or may put another
 *    idendifier which changes the calling convention of the function. If you
 *    really have to, ALWAYS declare it for the one platform you are dealing
 *    with:
 *
 *    Not ok:
 *        double exp(double a);
 *
 *    Ok:
 *        #ifdef SYMBOL_DEFINED_WEIRD_PLATFORM
 *        double exp(double);
 *        #endif
 */

#include <Python.h>
#include <math.h>

#include "config.h"
#include "numpy/npy_math.h"

/*
 *****************************************************************************
 **                     BASIC MATH FUNCTIONS                                **
 *****************************************************************************
 */

/* Original code by Konrad Hinsen.  */
#ifndef HAVE_EXPM1
static double expm1(double x)
{
    double u = exp(x);
    if (u == 1.0) {
        return x;
    } else if (u-1.0 == -1.0) {
        return -1;
    } else {
        return (u-1.0) * x/log(u);
    }
}
#endif

#ifndef HAVE_LOG1P
static double log1p(double x)
{
    double u = 1. + x;
    if (u == 1.0) {
        return x;
    } else {
        return log(u) * x / (u - 1);
    }
}
#endif

#ifndef HAVE_HYPOT
static double hypot(double x, double y)
{
    double yx;

    x = fabs(x);
    y = fabs(y);
    if (x < y) {
        double temp = x;
        x = y;
        y = temp;
    }
    if (x == 0.)
        return 0.;
    else {
        yx = y/x;
        return x*sqrt(1.+yx*yx);
    }
}
#endif

#ifndef HAVE_ACOSH
static double acosh(double x)
{
    return 2*log(sqrt((x+1.0)/2)+sqrt((x-1.0)/2));
}
#endif

#ifndef HAVE_ASINH
static double asinh(double xx)
{
    double x, d;
    int sign;
    if (xx < 0.0) {
        sign = -1;
        x = -xx;
    }
    else {
        sign = 1;
        x = xx;
    }
    if (x > 1e8) {
        d = x;
    } else {
        d = sqrt(x*x + 1);
    }
    return sign*log1p(x*(1.0 + x/(d+1)));
}
#endif

#ifndef HAVE_ATANH
static double atanh(double x)
{
    if (x > 0) {
        return -0.5*log1p(-2.0*x/(1.0 + x));
    }
    else {
        return 0.5*log1p(2.0*x/(1.0 - x));
    }
}
#endif

#ifndef HAVE_RINT
static double rint(double x)
{
    double y, r;

    y = floor(x);
    r = x - y;

    if (r > 0.5) goto rndup;

    /* Round to nearest even */
    if (r==0.5) {
        r = y - 2.0*floor(0.5*y);
        if (r==1.0) {
        rndup:
            y+=1.0;
        }
    }
    return y;
}
#endif

#ifndef HAVE_TRUNC
static double trunc(double x)
{
    return x < 0 ? ceil(x) : floor(x);
}
#endif

#ifndef HAVE_EXP2
#define LOG2 0.69314718055994530943
static double exp2(double x)
{
    return exp(LOG2*x);
}
#undef LOG2
#endif

#ifndef HAVE_LOG2
#define INVLOG2 1.4426950408889634074
static double log2(double x)
{
    return INVLOG2*log(x);
}
#undef INVLOG2
#endif

/*
 *****************************************************************************
 **                     IEEE 754 FPU HANDLING                               **
 *****************************************************************************
 */
#if !defined(HAVE_DECL_SIGNBIT)
#include "_signbit.c"

int _npy_signbit_f (float x)
{
    return _npy_signbit_d((double)x);
}

int _npy_signbit_ld (long double x)
{
    return _npy_signbit_d((double)x);
}
#endif

/*
 * if C99 extensions not available then define dummy functions that use the
 * double versions for
 *
 * sin, cos, tan
 * sinh, cosh, tanh,
 * fabs, floor, ceil, rint, trunc
 * sqrt, log10, log, exp, expm1
 * asin, acos, atan,
 * asinh, acosh, atanh
 *
 * hypot, atan2, pow, fmod, modf
 *
 * We assume the above are always available in their double versions.
 *
 * NOTE: some facilities may be available as macro only  instead of functions.
 * For simplicity, we define our own functions and undef the macros. We could
 * instead test for the macro, but I am lazy to do that for now.
 */

/**begin repeat
 * #type = npy_longdouble, float#
 * #TYPE = NPY_LONGDOUBLE, FLOAT#
 * #c = l,f#
 * #C = L,F#
 */

/**begin repeat1
 * #kind = sin,cos,tan,sinh,cosh,tanh,fabs,floor,ceil,rint,trunc,sqrt,log10,
 *         log,exp,expm1,asin,acos,atan,asinh,acosh,atanh,log1p,exp2,log2#
 * #KIND = SIN,COS,TAN,SINH,COSH,TANH,FABS,FLOOR,CEIL,RINT,TRUNC,SQRT,LOG10,
 *         LOG,EXP,EXPM1,ASIN,ACOS,ATAN,ASINH,ACOSH,ATANH,LOG1P,EXP2,LOG2#
 */

#ifdef @kind@@c@
#undef @kind@@c@
#endif
#ifndef HAVE_@KIND@@C@
static @type@ @kind@@c@(@type@ x)
{
    return (@type@) @kind@((double)x);
}
#endif

/**end repeat1**/

/**begin repeat1
 * #kind = atan2,hypot,pow,fmod#
 * #KIND = ATAN2,HYPOT,POW,FMOD#
 */
#ifdef @kind@@c@
#undef @kind@@c@
#endif
#ifndef HAVE_@KIND@@C@
static @type@ @kind@@c@(@type@ x, @type@ y)
{
    return (@type@) @kind@((double)x, (double) y);
}
#endif
/**end repeat1**/

#ifdef modf@c@
#undef modf@c@
#endif
#ifndef HAVE_MODF@C@
static @type@ modf@c@(@type@ x, @type@ *iptr)
{
    double niptr;
    double y = modf((double)x, &niptr);
    *iptr = (@type@) niptr;
    return (@type@) y;
}
#endif

/**end repeat**/

/* 
 * Useful constants in three precisions: 
 * XXX: those should really be in the header
 */

/**begin repeat
 * #c = f, ,l#
 * #C = F, ,L#
 */
#define NPY_E@c@        2.7182818284590452353602874713526625@C@ /* e */
#define NPY_LOG2E@c@    1.4426950408889634073599246810018921@C@ /* log_2 e */
#define NPY_LOG10E@c@   0.4342944819032518276511289189166051@C@ /* log_10 e */
#define NPY_LOGE2@c@    0.6931471805599453094172321214581766@C@ /* log_e 2 */
#define NPY_LOGE10@c@   2.3025850929940456840179914546843642@C@ /* log_e 10 */
#define NPY_PI@c@       3.1415926535897932384626433832795029@C@ /* pi */
#define NPY_PI_2@c@     1.5707963267948966192313216916397514@C@ /* pi/2 */
#define NPY_PI_4@c@     0.7853981633974483096156608458198757@C@ /* pi/4 */
#define NPY_1_PI@c@     0.3183098861837906715377675267450287@C@ /* 1/pi */
#define NPY_2_PI@c@     0.6366197723675813430755350534900574@C@ /* 2/pi */
/**end repeat**/

/*
 * Non standard functions
 */

/**begin repeat
 * #type = float, double, npy_longdouble#
 * #c = f, ,l#
 * #C = F, ,L#
 */

#define LOGE2    NPY_LOGE2@c@
#define LOG2E    NPY_LOG2E@c@
#define RAD2DEG  (180.0@c@/NPY_PI@c@)
#define DEG2RAD  (NPY_PI@c@/180.0@c@)

static @type@ rad2deg@c@(@type@ x) 
{
    return x*RAD2DEG;
}

static @type@ deg2rad@c@(@type@ x) 
{
    return x*DEG2RAD;
}

static @type@ log2_1p@c@(@type@ x)
{
    @type@ u = 1 + x;
    if (u == 1) {
        return LOG2E*x;
    } else {
        return npy_log2@c@(u) * x / (u - 1);
    }
}

static @type@ exp2_1m@c@(@type@ x)
{
    @type@ u = exp@c@(x);
    if (u == 1.0) {
        return LOGE2*x;
    } else if (u - 1 == -1) {
        return -LOGE2;
    } else {
        return (u - 1) * x/npy_log2@c@(u);
    }
}

static @type@ logaddexp@c@(@type@ x, @type@ y)
{
    const @type@ tmp = x - y;
    if (tmp > 0) {
        return x + npy_log1p@c@(npy_exp@c@(-tmp));
    }
    else {
        return y + npy_log1p@c@(npy_exp@c@(tmp));
    }
}

static @type@ logaddexp2@c@(@type@ x, @type@ y)
{
    const @type@ tmp = x - y;
    if (tmp > 0) {
        return x + log2_1p@c@(npy_exp2@c@(-tmp));
    }
    else {
        return y + log2_1p@c@(npy_exp2@c@(tmp));
    }
}

#define degrees@c@ rad2deg@c@
#define radians@c@ deg2rad@c@

#undef LOGE2
#undef LOG2E
#undef RAD2DEG
#undef DEG2RAD

/**end repeat**/

/*
 * Decorate all the functions: those are the public ones
 */

/**begin repeat
 * #type = npy_longdouble,double,float#
 * #c = l,,f#
 */
/**begin repeat1
 * #kind = sin,cos,tan,sinh,cosh,tanh,fabs,floor,ceil,rint,trunc,sqrt,log10,
 *         log,exp,expm1,asin,acos,atan,asinh,acosh,atanh,log1p,exp2,log2,
 *         rad2deg,deg2rad,exp2_1m#
 */

@type@ npy_@kind@@c@(@type@ x)
{
    return @kind@@c@(x);
}

/**end repeat1**/

/**begin repeat1
 * #kind = atan2,hypot,pow,fmod,logaddexp,logaddexp2#
 */
@type@ npy_@kind@@c@(@type@ x, @type@ y)
{
    return @kind@@c@(x, y);
}
/**end repeat1**/

@type@ npy_modf@c@(@type@ x, @type@ *iptr)
{
    return modf@c@(x, iptr);
}

/**end repeat**/
